Quantum heat transfer: A Born Oppenheimer method 
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We develop a Born-Oppenheimer type formalism for the description of quantum thermal transport along 
hybrid nanoscale objects. Our formalism is suitable for treating heat transfer in the off-resonant regime, where 
e.g., the relevant vibrational modes of the interlocated molecule are high relative to typical bath frequencies, 
and at low temperatures when tunneling effects dominate. A general expression for the thermal energy current 
is accomplished, in the form of a generalized Landauer formula. In the harmonic limit this expression reduces 
to the standard Landauer result for heat transfer, while in the presence of nonlinearities multiphonon tunneling 
effects are realized. 
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Introduction. — Thermal transport in molecular objects has 
recently become a topic of major interest, relevant for de- 
signing electronic and mechanical nanoscale devices U], and 
for resolving mechanisms and pathways of energy flow in 
biomolecules J2|. In modelling such systems we typically 
consider an impurity object, a subsystem, e.g., an alkane 
molecule J^, bridging two thermal reservoirs, representing 
solids or large residues in a protein, maintained each at a 
fixed temperature. Various treatments have been developed 
for simulating the thermal^conduction properties of such struc- 
tures, either classically fl], or in the quantum regime ||5l01. 
Among these treatments are the generalized Langevin equa- 
tion method Q [jj], t he Kinetic-Boltzmann theory Jgt], mode 
coupling the ory 111 fj . the non-equilibrium Green's function 
technique (HQj]], classical and mixed classical-quantum 
10, molecular dynamics simulations, and exact quantum 
simulations on simplified models 113 1. 

The master equation techni que at weak system-bath cou- 
pling is of particular interest |14fl, allowing to obtain sim- 
ple analytical results in interesting limits [15], guiding exper- 
imentalists and motivating theoreticians in developing more 
detailed treatments 1 1311 . In this approach the heat current is 
described by sequential incoherent emission and absorption 
processes, relaying on a resonance condition. Thus, a finite 
conductance exists only when the frequencies of the two ther- 
mal reservoirs match the subsystem characteristic frequency. 
Nevertheless, in many systems the characteristic frequencies 
of the impurity object are high relative to the cutoff frequen- 
cies of the reservoirs. For instance, consider an electronic spin 
surrounded by nuclear spins subjected to an external field, a 
molecule of high vibrational frequency coupled to solids with 
low Debye frequencies, or a high-frequency heat source in- 
side a protein with low frequency bonds as thermometers iflrjll . 
Developing a general formalism that can treat such scenarios, 
providing simple analytical results and bringing in physical 
insight, is of a great importance. 

Here we describe a new formalism for treating quantum 
thermal transport in such non-resonant systems, where sub- 
system's frequencies, relevant for thermal transfer, are above 
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FIG. 1: A Scheme of our setup (top), including a subsystem, e.g., a 
molecular chain, connecting two solids. The bottom panel exempli- 
fies the vibrational spectra p(ui) of the isolated solids and molecule. 



the reservoirs spectral window, or the baths temperatures are 
low, below the subsystem energy spacing 11711 . Using a Born- 
Oppenheimer (BO) type approximation we develop a compact 
expression for the heat current in the form of a generalized 
Landauer formula 111 811 . For harmonic systems we recover the 
elastic Landauer formula. When nonlinear interactions persist 
multiphonon processes contribute to the thermal current. 

Model. — Consider a small subsystem, representing e.g., 
a molecule, placed in between two thermal reservoirs (e.g., 
solids) maintained each at a fixed temperatures T v (v = L, R), 
see Fig. Q] The total Hamiltonian is given by 



H = H S + H L + Hf 



V L + V R , 



(1) 



where Hs is the Hamiltonian of the subsystem and H v stands 
for the v heat bath. Vl (Vr) couples separately the subsystem 
and the left (right) reservoir. The subsystem and the two reser- 
voirs' Hamiltonians assume diagonal forms, and we consider 
a bipartite form, V v = SB,y. Here S is a subsystem operator 
and B v is an operator in terms of the v bath degrees of free- 
dom. In particular, it is useful to study two extreme realiza- 
tions for the subsystem. In the first case nonlinear interactions 
are incorporated by adopting a two-level system (TLS) model 



1 1411 . S = a x , Hs = ^cr z . In the second model the subsystem 



is represented by a single harmonic mode linearly coupled to 
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the baths, S = b + tf, H s = Qtfb. Here & t (b) is the bosonic 
creation (annihilation) operator. 

Adiabatic evolution and the Bom-Oppenheimer 
approximation. — Consider the nonresonant case where 
the subsystem characteristic frequencies are high relative 
to the frequencies of the reservoirs fnll . This implies a 
timescale separation as the subsystem dynamics is fast, while 
the bath motion is slow. A BO type approximation can thus 
be employed following two consecutive steps: First, the fast 
variable is considered: We solve the subsystem eigenproblem 
fixing the reservoirs configuration, acquiring a set of potential 
energy surfaces which parametrically depend on B u . In the 
second step we assume that the baths dynamics evolves on 
the ground potential surface, and solve the vibrational heat 
transfer problem, form L to R. Next we follow this procedure 
using the generic Hamiltonian ((T). Beginning with the fast 
contribution, we diagonalize 

H g = H s + SB; (B = B L +B R ), (2) 

and acquire the potential surface W. For example, for a TLS 
subsystem we resolve \g(B v )) = cq |0) + c\ |1) as the ground 
state of H g ; co,i are the superposition coefficients, functions 
of B and e, with the eigenenergy W ~ — yJ{e/2) 2 + B 2 = 
-e/2-B 2 /e+B i /e 3 +0(B 6 /e 5 ). For an harmonic oscillator 
model we exactly obtain W = § — including the zero- 
point motion. We assume next that the total density matrix is 
initially factorized, 

p(0) = \g(B v ))p B (0)(g(B„)\, (3) 

where p B {0) = p L X p R ; p v = /Tr v [e~^] is 

the equilibrium-canonical distribution function of the v bath. 
Time evolution is dictated by the Liouville equation, 

p{t) = e - im p(0)e im 

« \g(B u ))e- iHBot p B (0)e iBBot (g(B^\, (4) 

where the second step is justified under the BO approximation 
with the effective Hamiltonian 

H BO = H l + H r + W. (5) 

Thus, the reduced density matrix of the reservoirs, 
PB(t) =Trsp(t), where the trace is performed over the sub- 
system degrees of freedom, evolves as 

p B (t) = e- lHBOt p B (0)e iHBot . (6) 

In the present scheme we thus propagate the bath coordinates 
along the subsystem potential energy surface W, and an ex- 
plicit study of the subsystem motion is not required, unlike 
the typical situation in other approaches flSQ- We iden- 
tify the operator W as an interaction term directly connecting 
the two reservoirs. Note that in the original model, Eq. (fTJ, V 
is linear in B, additive in the L and R coordinates. In contrast, 
under the BO approximation we obtain a potential energy sur- 
face W which is often nonlinear in B, mixing the left and 
right reservoirs' coordinates in a nontrivial way. 



Heat current. — The heat current operator, between the two 
reservoirs, can be defined as lfl9ll 

J L =i[H L ,W]/2. (7) 

For example, for a harmonic subsystem we recover Jl = 

— ^(BPl + PlB), while for a TLS subsystem « 
~i- e (BP L + P L B) when B/e « 1; P L = i[H L ,B L }. The 
current operator in both cases is identical in the first order of 
B/e. Generally, the expectation value of the current is 

J L (t) = Tr[J LPB (t)} = Tr[e^ BOt J L e-^ BOt ps(0)], (8) 

where the left expression is written in the Schoredinger pic- 
ture; the second is in the Heisenberg representation. The trace 
is performed over the two baths degrees of freedom. 

First Order Current — When system-baths couplings, ab- 
sorbed in W, are weak, the time evolution operator can be 
approximated by the first order term 

e -iH BO t = e -i(H L +H R )t (l-i j W{r)dT^j , (9) 

and the current (H)) reduces to 

J L {t) = -i[ Tr{[J L (r),^] Pi3 (0)}dT, (10) 

where W(t) and Jl(t) are interaction picture operators, 
A(t) = e iH B tj^ e -iH B t with Hb _ h l + H r . We are mostly 
interested in steady state quantities, J = J^it — > oo), if the 
limit exists. This expression can be further customized by us- 
ing a diagonal form for the reservoirs, e.g., for the L bath we 
write, Hl = Ek \k) (k\, and by expanding the potential 
surface in the left bath (L) and right (R) bath operators, func- 
tions of B v , 

w = Y. La ® Rb = Y.Y.Y. L tmK,s \ k p) ( ms \ ■ (U) 

a,b a.b k.m p,s 

\k) and |m) are the many body states of the left reservoir 
with energies Ek and E m ; \p) and |s) are the many body 
states of the right reservoir with energies E p and E s . The in- 
teraction W sums (nonseparable) contributions from the two 
reservoirs, a and b are integers. For example, for the har- 
monic subsystem with bilinear coupling W = —B 2 /fl = 

— {B\ + B 2 R + 2BlB b ) /fl. It can be shown that terms con- 
taining either L or R operators do not add to the current, 
as only mixed terms account. Therefore, in the case of a 
harmonic subsystem a single term contributes to (fTTl) with 
L 1 = iB L yj2/tt and R 1 = iB R y / 2/U. Back to (QjB, em- 
ploying (fTTT i. we accede to the steady state heat current 

J = z z X] X] \ L km\ 2 \ Rb ps\ 

^ ^ a,b k,m,p,s 

x5(E km + E ps )e-^ E "-^ E - 

— 2?I " \ ^ (\ T+ a I 2 I R- b | 2 I T~ a I 2 I R+ b \ 2 \ 

~ 2 Z ' ' ' ' ^' fcm ' \ U ps\ ~ l-^feml \ n ps | J 
a.b k,m,p,s 

x \E km \S{\E krn \ - \E ps \)e-^ E *-^ E ?, (12) 
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where E km = E k - E m and, e.g., Z L = J2 k e 0lE " is tne 
L bath partition function. Here L^m C^fem) denotes matrix 
elements when E k > E m (E k < E m ). We identify next the 
Fermi-like golden rule excitation (+) and relaxation (— ) rates, 
e.g., at the L contact, by 



fc±» = 27r]T \Lt a J 2 S(E k - E mT cj)- 

k,m 



(13) 



satisfying detailed balance, k^ a (oj) = k L a (ui)e ^ LLU . We can 
therefore reduce Eq. ( fT2b into the compact form 



1 f°° 

Z7T a,b J " 



-Pru> 



)(14) 



where the sum over a and b is determined given a particu- 
lar W. This is the main result of our paper. We refer to 
this expression as the "generalized Landauer formula" 11811 . 
as the net heat current is given by the difference between left- 
moving and right-moving excitations. Nevertheless, our for- 
mula can incorporate anharmonic interactions, absorbed in the 
rates k^ a (u>) unlike the original treatment 1 1811 . We empha- 
size the broad status of Eq. (TT4-b . It was derived without spec- 
ifying the subsystem Hamiltonian or the system-bath interac- 
tion form, both contained in W. It is valid as long as (i) there 
exists a timescale separation between the subsystem motion 
(fast) and the reservoirs dynamics (slow), and (ii) system-bath 
interaction is weak, see Eq. (O. In what follows we apply 
Eq. ( TPfl i on some models of particular interest: a fully har- 
monic model, a nonlinear model with strong system-bath in- 
teractions, and utilizing a spin subsystem, representing a non- 
linear impurity. 

Harmonic model. — We consider first a harmonic model, 
H = H l +H r + Vl + V r + H s , with 

S=(&t + fc), /;, ^ A,, , ;/;,, , • /,:; ,;. (15) 

and show that Eq. ( fT4l i reduces to the standard elastic limit 
J3, O- Here the subsystem comprises a single mode of 
frequency Q, i>£ . (b v j) are the creation (annihilation) op- 
erators of the mode j in the v bath, and b are the re- 
spective subsystem operators. A„j are system-bath interac- 
tion energies, S is a subsystem operator. The expectation 
value of the current can be calculated either by following 
Eq. ( TTOb in the long time limit, or by directly applying Eq. 
(O, as we do next. In the occupation number represen- 
tation the many body states of the L reservoir are | to) = 
|toi, m2...m;...mjv) with to; excitations for the I mode. Since 
W = —B 2 /tt, the relevant matrix elements in (fTTb are 

\ L Li\ = \[h E; x lj (Vm + is kl , mi +i + ^a-i). 

An analogous expression exists for Ri s . We thus identify 
the excitation and relaxation rates in Eq. (IT3b by fc+ 1 (w) = 



^r v (u})n u (uj) and k y x {uj) = 

tively, where n v {uS) = [e w ' T,y 
distribution function and T v (ui) 



n r„(w) [n v {u) 



1 



1], respec- 



is the Bose-Einstein 
Using these rates the current (TBI) reduces to 



J 



T L (u>)T R (u>) 

n 2 



[nL{w) — n R (u>)} uduj. (16) 



This is the Landauer's formula for heat conduction [18] in the 
BO limit; assuming the system frequency is above the baths 
spectral window, further utilizing the weak-coupling approx- 
imation [Eq. (|9|. Beyond this limit, the heat current for the 
harmonic model (TT3T > is exactly given by 

^ T{uj)[n L (uj) - n R (uj)]ujduj, (17) 



.7 



with 



the 

fined above Eq. (16 



transmission 

fzh. 



coefficient T(u>) = 
The rate r„ has been de- 



for convenience we discard the direct 
reference to frequency. In the weak coupling limit, T u < 0, 
the transmission coefficient is sharply peaked around f2. In 
the nonresonant case f2 3> lo c , where u> c is the reservoirs 
cutoff frequency, T(cj) ~ ^M^rM and Eq ^ is 
recovered. In the opposite limit, when the baths spectral 
window overlap with the molecular vibrations, us c ^> fi, 
Eq. ( TPTI i reduces into a resonant energy transfer expression, 
J = ^ r^r ["l(^) - n R (Q,)]. Here T v is calculated at the 
(local oscillator) frequency f2. This expression describes a 
hopping motion, with energy flowing sequentially from the L 
bath into the subsystem, then into the R contact. This process 
is dictated by the subsystem energetic window, yielding 
J oc f2. In contrast, Eq. ( fT6b accounts for a coherent, deep 
tunneling energy transfer mechanism, and the current decays 
with the energetic barrier, J oc I/O 2 . 

Anharmonic models. — We generalize next the harmonic re- 
sult by modifying the model (fT5l ). adopting an exponentially 
repulsive interaction, 

e -E,-M&t iJ+ &„, f ) 



(18) 



appropriate for the relevant nonresonant case 12011 . As before, 
diagonalizing H g = Hs + Vl + Vr we obtain the poten- 
tial surface W = -B 2 /Q; B = B L + Br, controlled by 
the bipartite term 2BlB r /Q. For simplicity, we assume an 
(identical) Einstein-type model for the reservoirs spectra, rep- 
resented by a single frequency lj r ■ Under this assumption the 
relevant excitation/relaxation rates [Eq. ( fOb l are given by 



9 \ 21 71 

E 2 4EaAyM^) + ir 

\l-e 



kt\w) 



fi 



x n v {u)B) s (5(t^ - (2s - 1)wb)- 



(19) 



Assuming A is small, we enclose only single-phonon and two- 
phonon contributions in (fl4] i. yielding the heat current 

J = ^{wflA|A|(ni - n R ) 
, (2^ B )A 4 L A 4 fl 



[n 2 L (n R + I) 2 - n 2 R (n L + I) 2 ] }.(20) 
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The Bose-Einstein functions are evaluated at the frequency 
u>B- This expression presents a generalization to the harmonic 
result ( fT6] l, accommodating multiphonon processes; the sec- 
ond term describes tunneling of a two-phonon combination. 

The starting point of our next anharmonic model is again 
Eq. ( TTBT l. utilizing a two-level subsystem, H$ = ^<J Z , S = 
a x , representing a nonlinear impurity bilinearly coupled to 
bath phonons. In this case we resolve W = — V B 2 + eJ 4, 
thus the first order current is identical to the harmonic result 
(fl6l l. with il replaced by e. Incorporating the next term in the 
expansion, — B 2 /e + B 4 /e 3 , we get 

1 f°° 

J=— uoduj^- 13 ^ -e~ PR ")[kl 1 (uj)k- R 1 {uj) + 
2 -71 " Jo 

ki 2 (uj)k R 2 {u:) + fc^H^'M + fcz'Mfc^M]^) 

The first term in the square brackets describes a single phonon 
(harmonic) process, proportional to 1/e 2 . The other terms 
collect contributions from multiphonon processes. For exam- 
ple, the second element accounts for the absorption of two 
phonons in the left bath, followed by an emission of these 
phonons at the other end, with, e.g., 

fc^ 2 (w) oc 4^ ^(1 + + ni')8{uj - u>i - uj v ) 
i.v 

+ \^2{l + ni )n v 5{u -ui+u v ). (22) 

U' 

The last two contributions in d2"TT i convene three-phonon pro- 
cesses, where, e.g., a single mode from the left bath decays 
into three excitations at the right side. Fig. [2] presents the 
frequency components of the heat current (the integrand of 
Eq. (ETJ), where for simplicity we assume a spectral density 
S(lo) = J2j j3( w — w j) P ea ked around a specific bath fre- 
quency lob = 2 with a hard cutoff at oj = 3, see panel (a). 
We identify three contributions to the current: A dominant, 
single-phonon element at u> ~ uj b , a weaker two-phonon con- 
tribution, and a rudimentary three-phonon current, see panel 
(b). In the presence of a spatial asymmetry these high or- 
der terms are responsible for the thermal rectification effect 
1 2 lfl . Note that in the resonant regime, when a hopping mech- 
anism dominates, the heat current across harmonic junctions 
is higher than its anharmonic counterpart II 1411 due to a satura- 
tion effect. In contrast, in the nonresonant case anharmonicity 
enhances the thermal current due to the participation of mul- 
tiphonon processes. Similar observations were obtained in a 



study of classical heat flow in molecular junctions [22]. 

Summary. — We have presented here a generally applica- 
ble Born-Oppenheimer type formalism for describing thermal 
energy transfer in the off -resonant case, where an impurity ob- 
ject has a characteristic frequency above the (populated) bath 
modes. In this limit energy propagates across the structure in 
a tunneling-like motion, keeping the subsystem population in- 
tact. In the weak coupling limit we derived a compact expres- 
sion for the thermal current, bearing the structure of a general- 
ized Landauer relation, yet incorporating multiphonon effects. 




FIG. 2: Frequency components of the heat current J(u>), for a spin 
subsystem bilinearly coupled to heat baths. Multiphonon processes 
are observed, see also panel (b). The parameters e = 6, Tl = 1, 
Tr = 0.5 were used with the bath spectral function S(lo), depicted 
in panel (b), identical at the two ends. 



In the harmonic limit our formula reduces to known results. 
We have also applied it onto nonlinear models: Incorporating 
molecular anharmonicity or assuming short range interactions 
we reach simple analytic expressions for the heat current, re- 
flecting the underling transport mechanism. The new method 
described here is complementary to kinetic approaches that 
are typically valid in the resonant case. By incorporating 
quantum effects and nonlinearity it opens new ways for de- 
scribing vibrational or electronic energy dynamics in organic 
materials 1231. biomolecules J2l and superconductors |24]. 
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